function expandFromPoint(inputFileName,sliceNumber,seed_diameter,outputFileName)

%  tb3Jc2Ee2WW8
%input:
%inputFileName - original image file name
%sliceNumber - number of slice to look for the seed
%seed_diameter - seed diameter
%outputFileName - filtered image save path

if (nargin~=4)
    error('wrong number of aruments');  
end

close all;
seed_diameter = floor(seed_diameter/2);

%% Load and display the file
struct = load_untouch_nii(inputFileName);
srcImg = struct.img;
resImg = zeros(size(srcImg)); % output image

%% pool the picture, show on screen and wait for user input

h=figure(1);
imshow(srcImg(:,:,sliceNumber),[]);
[userY, userX] = ginput(1); % read one user mouse input
close(h);
userX = round(userX);
userY = round(userY);

userPixelVal =srcImg(userY,userX,sliceNumber);
[Xn, Yn, Zn] = size(srcImg);



queue =  zeros(Yn*Xn*Zn,3);
regionValues=zeros(Xn*Yn*Zn*3,1);
cursorTreshold=1;
cursor=1;
for i = -seed_diameter:seed_diameter
    for j = -seed_diameter:seed_diameter
        for k = sliceNumber-seed_diameter: sliceNumber +seed_diameter            
            x = userX + i;
            y = userY +j;
            z = k;            
            if isPointInsideImage(x,y,z, srcImg) %check if the point is inside image dimensions.
                queue(cursor,:) = [x, y, z];         
                regionValues(cursorTreshold) = srcImg(y,x,z);
                cursorTreshold=cursorTreshold+1;
                cursor =  cursor+1;                
            end
        end
    end
end


treshold = abs( (max(regionValues(1:cursorTreshold-1))) - (min(regionValues(1:cursorTreshold-1)))) ;


cursorNow=1;
cursorNext=cursor; % next empty row
while( cursorNow < cursorNext )
   % get next voxel value from Queue
    x = queue(cursorNow,1);
    y = queue(cursorNow,2);
    z = queue(cursorNow,3);
    queue(cursorNow,:) = [0 0 0];    
    currVal = srcImg(y,x,z);
    resImg(y,x,z) = 1; % add "testing voxel" to the result image.
     
    % check neighbours from current voxel.    
    for i = -1:2:1
        for j = -1:2:1
            for k = -1:2:1                
                    x=x+i;y=y+j;z=z+k;
                    if (isPointInsideImage(x,y,z,srcImg))                        
                        %isInRange = srcImg(y,x,z) <= (userPixelVal + treshold) && srcImg(y,x,z) >= (userPixelVal-treshold);
          
                        isInRange =  srcImg(y,x,z) >= (userPixelVal-treshold) && srcImg(y,x,z) <= (userPixelVal + treshold) ;
                        if (isInRange==1 &&  resImg(y,x,z) ==0 )
                            regionValues(cursorTreshold) = srcImg(y,x,z);
                            cursorTreshold=cursorTreshold+1;
                            queue(cursorNext,:) = [x,y,z]; % add neighbour to testing pixel.
                            cursorNext=cursorNext+1;
                            resImg(y,x,z) = 1; % mark the voxel as processed.
                           % treshold = std(regionValues(1:cursorTreshold-1,:));
                            %imshow(resImg(:,:,sliceNumber),[]);
                        end
                    end                
            end
        end
    end    
 % treshold = std(regionValues(1:cursorTreshold-1,:));
    %treshold = double( (max(regionValues(:,:))) - (min(regionValues(:,:)))) * 0.05;
    treshold = abs( (max(regionValues(1:cursorTreshold-1))) - (min(regionValues(1:cursorTreshold-1)))) ;
    cursorNow = cursorNow+1;
    display(cursorNow);
end %end while

%Finish
%sub2ind
%bwmorph
%region

struct.img = resImg;
save_untouch_nii(struct,outputFileName);

%%
figure(1);
imshow(srcImg(:,:,sliceNumber),[]);
figure(2);
imshow(resImg(:,:,sliceNumber),[]);
disp('expandFromPoint --> FINISHED');
end


function  res = isPointInsideImage(x,y,z,img)
% checks if current point is inside image dimensions
    [Xn Yn Zn ] = size(img);
    res = ( x>0 && x<=Xn) && (y>0 && y<=Yn) && (z>0 && z<=Zn); 
end

